In [1]:
using ValidatedNumerics

In [2]:
logistic(r, x) = -r*((x-0.5)^2 - 0.25) #r*x*(1-x)


Out[2]:
logistic (generic function with 1 method)

In [3]:
X = @interval(-0.0001, 1.0001)


Out[3]:
[-0.0001, 1.0001000000000002]

In [4]:
r = @interval(3.1)
newton(x -> logistic(r,x)-x, X)


Out[4]:
2-element Array{(Interval{Float64},Symbol),1}:
 ([-6.968514798756479e-17, 1.3517743393731552e-16],:unique)
 ([0.6774193548387095, 0.6774193548387097],:unique)        

In [5]:
r = @interval(3.1)
newton(x -> logistic(r, logistic(r,x))-x, @interval(-1, 2))


Out[5]:
4-element Array{(Interval{Float64},Symbol),1}:
 ([-6.233702290330726e-17, 1.1931689491945523e-16],:unique)
 ([0.5580141252026949, 0.5580141252026971],:unique)        
 ([0.6774193548387076, 0.6774193548387121],:unique)        
 ([0.7645665199585938, 0.7645665199585949],:unique)        

In [3]:
function iterate_logistic(n)
    ex = :(x)
    for i in 1:n
        ex = :(logistic(r,$ex))
    end
    ex = :($(ex) - x)
    ex
end

for n in [1, 2, 4, 5, 8:16, 32]
    name = symbol("logistic_iterated_$n")
    ex = iterate_logistic(n)
    @eval $name(r,x) = $ex
end

In [5]:
iterate_logistic(4)


Out[5]:
:(logistic(r,logistic(r,logistic(r,logistic(r,x)))) - x)

In [14]:
logistic_iterated_16(3.1, 0)


Out[14]:
-0.0

In [23]:
r = @interval(3.1)
newton(x -> logistic_iterated_2(r, x), X)


Out[23]:
4-element Array{(Interval{Float64},Symbol),1}:
 ([-9.947150462130687e-17, 1.1038772564320529e-16],:unique)
 ([0.5580141252026949, 0.5580141252026971],:unique)        
 ([0.6774193548387076, 0.6774193548387116],:unique)        
 ([0.7645665199585937, 0.7645665199585948],:unique)        

In [22]:
r = @interval(3.5)
newton(x -> logistic_iterated_4(r, x), X)


Out[22]:
8-element Array{(Interval{Float64},Symbol),1}:
 ([-9.361862579270038e-17, 1.0515086331210236e-16],:unique)
 ([0.3828196830173237, 0.3828196830173243],:unique)        
 ([0.4285714285714273, 0.4285714285714295],:unique)        
 ([0.5008842103072174, 0.5008842103072186],:unique)        
 ([0.7142857142857141, 0.7142857142857144],:unique)        
 ([0.8269407065914381, 0.8269407065914388],:unique)        
 ([0.8571428571428567, 0.8571428571428578],:unique)        
 ([0.874997263602464, 0.8749972636024642],:unique)         

In [21]:
r = @interval(3.55)
newton(x -> logistic_iterated_8(r, x), X)


Out[21]:
16-element Array{(Interval{Float64},Symbol),1}:
 ([-6.450867901675527e-17, 1.5719872312603065e-16],:unique)
 ([0.3548004479999936, 0.35480044799999727],:unique)       
 ([0.36112003982837154, 0.36112003982838603],:unique)      
 ([0.37032556106597025, 0.37032556106597875],:unique)      
 ([0.41803814056962735, 0.41803814056962996],:unique)      
 ([0.5060305096360181, 0.5060305096360388],:unique)        
 ([0.5261830680927085, 0.5261830680927377],:unique)        
 ([0.5404748339895925, 0.5404748339896022],:unique)        
 ([0.7183098591549292, 0.71830985915493],:unique)          
 ([0.8126556698514403, 0.8126556698514448],:unique)        
 ([0.8190288661526807, 0.8190288661526939],:unique)        
 ([0.8278051165993622, 0.8278051165993711],:unique)        
 ([0.8636520002754412, 0.8636520002754425],:unique)        
 ([0.8816843467379748, 0.8816843467379784],:unique)        
 ([0.8850662866556411, 0.8850662866556468],:unique)        
 ([0.8873708969850299, 0.8873708969850314],:unique)        

In [20]:
r = @interval(3.55)
krawczyk(x -> logistic_iterated_8(r, x), X)


Out[20]:
16-element Array{(Interval{Float64},Symbol),1}:
 ([-9.468711597088143e-17, 1.0093800912363333e-16],:unique)
 ([0.35480044799999333, 0.3548004479999971],:unique)       
 ([0.3611200398283715, 0.3611200398283861],:unique)        
 ([0.3703255610659702, 0.3703255610659793],:unique)        
 ([0.41803814056962735, 0.41803814056963],:unique)         
 ([0.5060305096360179, 0.5060305096360391],:unique)        
 ([0.5261830680927094, 0.5261830680927403],:unique)        
 ([0.5404748339895923, 0.5404748339896025],:unique)        
 ([0.7183098591549291, 0.7183098591549301],:unique)        
 ([0.8126556698514402, 0.8126556698514449],:unique)        
 ([0.8190288661526809, 0.8190288661526951],:unique)        
 ([0.8278051165993631, 0.8278051165993718],:unique)        
 ([0.8636520002754411, 0.8636520002754426],:unique)        
 ([0.8816843467379746, 0.8816843467379786],:unique)        
 ([0.8850662866556408, 0.8850662866556466],:unique)        
 ([0.8873708969850298, 0.8873708969850315],:unique)        

In [19]:
r = @interval(4)
newton(x -> logistic_iterated_8(r, x), X)


Out[19]:
256-element Array{(Interval{Float64},Symbol),1}:
 ([-6.685942182571293e-17, 9.736289924187537e-17],:unique) 
 ([0.00014942107845312833, 0.0001494210784532221],:unique) 
 ([0.00015177401106415282, 0.00015177401106424558],:unique)
 ([0.0005975950071779156, 0.0005975950071780078],:unique)  
 ([0.0006070039028550085, 0.0006070039028551082],:unique)  
 ([0.0013442539196471853, 0.001344253919647281],:unique)   
 ([0.0013654133071059132, 0.0013654133071060027],:unique)  
 ([0.0023889515495412795, 0.0023889515495413996],:unique)  
 ([0.00242654179646774, 0.002426541796467861],:unique)     
 ([0.0037310634974740595, 0.0037310634974741606],:unique)  
 ([0.0037897451640320594, 0.003789745164032159],:unique)   
 ([0.005369787604186948, 0.005369787604187054],:unique)    
 ([0.00545419581442697, 0.005454195814427069],:unique)     
 ⋮                                                         
 ([0.9954156265457801, 0.9954156265457802],:unique)        
 ([0.9954865939027027, 0.9954865939027029],:unique)        
 ([0.9969295684476368, 0.9969295684476369],:unique)        
 ([0.9969771232924401, 0.9969771232924403],:unique)        
 ([0.9981418263742147, 0.9981418263742148],:unique)        
 ([0.9981706172512619, 0.998170617251262],:unique)         
 ([0.999051664368522, 0.9990516643685221],:unique)         
 ([0.9990663624465502, 0.9990663624465503],:unique)        
 ([0.9996585300715114, 0.9996585300715115],:unique)        
 ([0.9996638235054527, 0.9996638235054528],:unique)        
 ([0.9999620550574152, 0.9999620550574153],:unique)        
 ([0.9999626433348662, 0.9999626433348664],:unique)        

In [24]:
r = @interval(4)
newton(x -> logistic_iterated_16(r, x), X)


Out[24]:
220598-element Array{(Interval{Float64},Symbol),1}:
 ([-0.0001, -9.999906849116088e-5],:unknown)              
 ([-9.999906849116087e-5, -9.999813698232174e-5],:unknown)
 ([-9.999813698232173e-5, -9.99972054734826e-5],:unknown) 
 ([-9.999720547348259e-5, -9.999627396464349e-5],:unknown)
 ([-9.999627396464347e-5, -9.999534245580435e-5],:unknown)
 ([-9.999534245580434e-5, -9.999441094696521e-5],:unknown)
 ([-9.99944109469652e-5, -9.999347943812607e-5],:unknown) 
 ([-9.999347943812606e-5, -9.999254792928696e-5],:unknown)
 ([-9.999254792928695e-5, -9.999161642044782e-5],:unknown)
 ([-9.99916164204478e-5, -9.999068491160868e-5],:unknown) 
 ([-9.999068491160867e-5, -9.998975340276954e-5],:unknown)
 ([-9.998975340276953e-5, -9.998882189393043e-5],:unknown)
 ([-9.998882189393042e-5, -9.998789038509129e-5],:unknown)
 ⋮                                                        
 ([1.0000999888218949, 1.0000999897534038],:unknown)      
 ([1.000099989753404, 1.0000999906849124],:unknown)       
 ([1.0000999906849126, 1.000099991616421],:unknown)       
 ([1.0000999916164213, 1.0000999925479297],:unknown)      
 ([1.00009999254793, 1.0000999934794388],:unknown)        
 ([1.000099993479439, 1.0000999944109474],:unknown)       
 ([1.0000999944109477, 1.0000999953424565],:unknown)      
 ([1.0000999953424567, 1.0000999962739652],:unknown)      
 ([1.0000999962739654, 1.0000999972054743],:unknown)      
 ([1.0000999972054745, 1.000099998136983],:unknown)       
 ([1.0000999981369831, 1.0000999990684916],:unknown)      
 ([1.0000999990684918, 1.0001000000000002],:unknown)      

In [26]:
roots = ans;
filter(x -> x[2]==:unique, roots)


Out[26]:
65069-element Array{(Interval{Float64},Symbol),1}:
 ([1.3273284430570427e-5, 1.327328443066284e-5],:unique) 
 ([1.7795734016136777e-5, 1.779573401628048e-5],:unique) 
 ([1.8612681896730634e-5, 1.8612681896823096e-5],:unique)
 ([2.0305136135843612e-5, 2.030513613593623e-5],:unique) 
 ([2.1620561528009657e-5, 2.1620561528102496e-5],:unique)
 ([2.3440452157732228e-5, 2.3440452157831687e-5],:unique)
 ([2.344188288191098e-5, 2.344188288201578e-5],:unique)  
 ([2.4379458657892453e-5, 2.4379458657984908e-5],:unique)
 ([2.4855140261344045e-5, 2.4855140261437744e-5],:unique)
 ([2.5820290877229384e-5, 2.5820290877321853e-5],:unique)
 ([2.6309759880782993e-5, 2.630975988088592e-5],:unique) 
 ([2.680218876442067e-5, 2.6802188764513737e-5],:unique) 
 ([2.730081891125118e-5, 2.7300818911420588e-5],:unique) 
 ⋮                                                       
 ([0.9999734437821953, 0.9999734437821955],:unique)      
 ([0.9999734454029944, 0.9999734454029946],:unique)      
 ([0.9999739371398805, 0.9999739371398807],:unique)      
 ([0.9999744227202014, 0.9999744227202015],:unique)      
 ([0.9999744242812537, 0.9999744242812538],:unique)      
 ([0.9999749052955099, 0.99997490529551],:unique)        
 ([0.9999753847774436, 0.9999753847774437],:unique)      
 ([0.9999758581322514, 0.9999758581322516],:unique)      
 ([0.9999763254466042, 0.9999763254466043],:unique)      
 ([0.9999772506234746, 0.9999772506234748],:unique)      
 ([0.9999879210995385, 0.9999879210995387],:unique)      
 ([0.9999879218367506, 0.9999879218367508],:unique)      

In [34]:
logistic_iterated_16(3.1, 0)


Out[34]:
0.0

In [29]:
newton(x->logistic_iterated_5(r, x), X)


Out[29]:
32-element Array{(Interval{Float64},Symbol),1}:
 ([-5.660020044625717e-17, 1.3150105603012344e-16],:unique)
 ([0.00903565136864659, 0.009035651368646704],:unique)     
 ([0.010235029373752695, 0.010235029373752817],:unique)    
 ([0.0358160334919636, 0.03581603349196378],:unique)       
 ([0.040521094189884616, 0.0405210941898848],:unique)      
 ([0.07937323358440937, 0.07937323358440947],:unique)      
 ([0.08961827939636176, 0.08961827939636192],:unique)      
 ([0.1381329809474648, 0.13813298094746498],:unique)       
 ([0.15551654046215668, 0.15551654046215696],:unique)      
 ([0.20997154521440084, 0.20997154521440098],:unique)      
 ([0.2355179948365187, 0.2355179948365189],:unique)        
 ([0.29229249349905667, 0.2922924934990569],:unique)       
 ([0.32634737357758975, 0.32634737357758997],:unique)      
 ⋮                                                         
 ([0.7201970757788171, 0.7201970757788173],:unique)        
 ([0.7499999999999999, 0.7500000000000001],:unique)        
 ([0.8060529912738313, 0.8060529912738316],:unique)        
 ([0.8274303669726424, 0.8274303669726426],:unique)        
 ([0.8793790613463953, 0.8793790613463955],:unique)        
 ([0.8930265473713936, 0.8930265473713939],:unique)        
 ([0.9371733080722909, 0.9371733080722912],:unique)        
 ([0.9444177243274616, 0.9444177243274618],:unique)        
 ([0.9770696282000243, 0.9770696282000245],:unique)        
 ([0.9797464868072486, 0.9797464868072487],:unique)        
 ([0.9974346616959475, 0.9974346616959476],:unique)        
 ([0.9977359612865422, 0.9977359612865424],:unique)        

In [30]:
roots[roots[2]==:unique]


BoundsError()
while loading In[30], in expression starting on line 1

 in getindex at array.jl:246

In [32]:
r = @interval(4)
newton(x -> logistic_iterated_10(r, x), X)


Out[32]:
1024-element Array{(Interval{Float64},Symbol),1}:
 ([-1.0220696221544258e-16, 8.283033843053199e-17],:unique)
 ([9.394002137712196e-6, 9.394002137833982e-6],:unique)    
 ([9.430769118680364e-6, 9.430769118773784e-6],:unique)    
 ([3.757565556196494e-5, 3.7575655562059124e-5],:unique)   
 ([3.772272071721069e-5, 3.772272071733585e-5],:unique)    
 ([8.454390131859121e-5, 8.454390131868317e-5],:unique)    
 ([8.487478753608691e-5, 8.487478753617993e-5],:unique)    
 ([0.00015029697452838232, 0.0001502969745285893],:unique) 
 ([0.0001508851908543323, 0.0001508851908544388],:unique)  
 ([0.00023483240445339927, 0.00023483240445349297],:unique)
 ([0.0002357514405563716, 0.00023575144055646552],:unique) 
 ([0.0003381470145895322, 0.00033814701458962634],:unique) 
 ([0.00033947033522630486, 0.0003394703352264008],:unique) 
 ⋮                                                         
 ([0.999714745464859, 0.9997147454648591],:unique)         
 ([0.9997158574619074, 0.9997158574619076],:unique)        
 ([0.9998090384816738, 0.999809038481674],:unique)         
 ([0.999809782923107, 0.9998097829231072],:unique)         
 ([0.9998844771639072, 0.9998844771639074],:unique)        
 ([0.9998849275276295, 0.9998849275276297],:unique)        
 ([0.9999410586657799, 0.9999410586657801],:unique)        
 ([0.9999412884518407, 0.9999412884518408],:unique)        
 ([0.9999787808528637, 0.9999787808528638],:unique)        
 ([0.9999788635779219, 0.9999788635779221],:unique)        
 ([0.9999976423021615, 0.9999976423021617],:unique)        
 ([0.99999765149395, 0.9999976514939501],:unique)          

In [33]:
length(filter(x->x[2]==:unique, ans))


Out[33]:
1024

In [40]:
r = @interval(4)
newton(x -> logistic_iterated_11(r, x), X, tolerance=0, maxlevel=40)


Out[40]:
2048-element Array{(Interval{Float64},Symbol),1}:
 ([-3.7739729411576187e-17, 7.443905684961557e-17],:unique)
 ([2.350798951126237e-6, 2.3507989511639162e-6],:unique)   
 ([2.3553948387492456e-6, 2.3553948388366314e-6],:unique)  
 ([9.403173699520289e-6, 9.403173699560357e-6],:unique)    
 ([9.421557163670454e-6, 9.421557163730942e-6],:unique)    
 ([2.1157057930306742e-5, 2.115705793034727e-5],:unique)   
 ([2.11984204002431e-5, 2.119842040030646e-5],:unique)     
 ([3.761234111944355e-5, 3.76123411194866e-5],:unique)     
 ([3.768587359185487e-5, 3.768587359191397e-5],:unique)    
 ([5.8768868534615475e-5, 5.8768868534701886e-5],:unique)  
 ([5.888376140063535e-5, 5.888376140069161e-5],:unique)    
 ([8.46264412368879e-5, 8.462644123696015e-5],:unique)     
 ([8.479188410901219e-5, 8.47918841090595e-5],:unique)     
 ⋮                                                         
 ([0.9999287509423916, 0.9999287509423918],:unique)        
 ([0.9999288899616146, 0.9999288899616147],:unique)        
 ([0.9999523039753886, 0.9999523039753887],:unique)        
 ([0.9999523970393035, 0.9999523970393037],:unique)        
 ([0.9999711466680797, 0.99997114666808],:unique)          
 ([0.9999712029667102, 0.9999712029667103],:unique)        
 ([0.9999852788429373, 0.9999852788429374],:unique)        
 ([0.9999853075669987, 0.9999853075669988],:unique)        
 ([0.9999947003668137, 0.9999947003668139],:unique)        
 ([0.9999947107075408, 0.9999947107075409],:unique)        
 ([0.9999994111509435, 0.9999994111509436],:unique)        
 ([0.9999994122999167, 0.9999994122999168],:unique)        

In [35]:
length(filter(x->x[2]==:unique, ans))


Out[35]:
2047

In [36]:
newton(x->logistic_iterated_11(r, x), @interval("[0.9999994111193403, 0.9999994112365287]"))


Out[36]:
1-element Array{(Interval{Float64},Symbol),1}:
 ([0.9999994111509435, 0.9999994111509436],:unique)

In [38]:
newton(x->logistic_iterated_11(r, x), @biginterval("[0.9999994111193403, 0.9999994112365287]"))


Out[38]:
1-element Array{(Interval{BigFloat},Symbol),1}:
 ([9.999994111509435535683237673405515492356983378423455810550216538164422108648822e-01, 9.999994111509435535683237673405515492356983378423455810550216538164422108648995e-01]₂₅₆,:unique)

In [39]:
diam(ans[1][1])


Out[39]:
1.727233711018888925077270372560079914223200072887256277004740694033718360632485e-77 with 256 bits of precision

Exact periodic points for $r=4$

When $r=4$, it is possible to calculate the periodic points exactly; see e.g. Sanders+Benet, Scipy 2014 proceedings. First we calculate the periodic points of the tent map using BigFloats:


In [25]:
function tent_map_periodic_points(n)
    
    evens = big(1.0)/(big(2)^n-1) * collect(0:2:2^n-1)
    odds = big(1.0)/(big(2)^n+1) * collect(2:2:2^n)
    
    points = zeros(BigFloat, 2^n)
    points[1:2:2^n] = evens
    points[2:2:2^n] = odds
    
    points
    
end


Out[25]:
tent_map_periodic_points (generic function with 1 method)

In [26]:
tent_map_periodic_points(4)


Out[26]:
16-element Array{BigFloat,1}:
 0e+00                                                                               
 1.17647058823529411764705882352941176470588235294117647058823529411764705882353e-01 
 1.333333333333333333333333333333333333333333333333333333333333333333333333333343e-01
 2.35294117647058823529411764705882352941176470588235294117647058823529411764706e-01 
 2.666666666666666666666666666666666666666666666666666666666666666666666666666687e-01
 3.529411764705882352941176470588235294117647058823529411764705882352941176470601e-01
 4.000000000000000000000000000000000000000000000000000000000000000000000000000052e-01
 4.70588235294117647058823529411764705882352941176470588235294117647058823529412e-01 
 5.333333333333333333333333333333333333333333333333333333333333333333333333333374e-01
 5.882352941176470588235294117647058823529411764705882352941176470588235294117683e-01
 6.666666666666666666666666666666666666666666666666666666666666666666666666666695e-01
 7.058823529411764705882352941176470588235294117647058823529411764705882352941202e-01
 8.000000000000000000000000000000000000000000000000000000000000000000000000000104e-01
 8.235294117647058823529411764705882352941176470588235294117647058823529411764721e-01
 9.333333333333333333333333333333333333333333333333333333333333333333333333333425e-01
 9.41176470588235294117647058823529411764705882352941176470588235294117647058824e-01 

These points $y_n$ are mapped onto periodic points of the logistic map by $x := h(y) := \sin^2(\frac{\pi}{2} y)$


In [27]:
h(y) = sin(π*y/2)^2


Out[27]:
h (generic function with 1 method)

In [28]:
logistic_exact_period_4_points = map(h, tent_map_periodic_points(4))


Out[28]:
16-element Array{BigFloat,1}:
 0e+00                                                                               
 3.376388529782209771344205408921830686870611102744153587582499406973197671517772e-02
 4.322727117869955224893621400734141102959477031126272746950010596744255894868408e-02
 1.304955413896704420377328450636759471200411901459511384214268446520135907271199e-01
 1.654346968205708930868633346566097632002083905201021615912733238560166991446783e-01
 2.771308221117308663017712253102565722861611132052862534998310635742086249501397e-01
 3.454915028125262879488532914085904705699227050485592844661378443236848842952812e-01
 4.538658202683490023801744464227467598184913581220883334042287350394312129893768e-01
 5.522642316338267356999170774012490595403279347372965569966377383562423253465847e-01
 6.368314950360414317695389677184067158124263347132527319387145577056778723541536e-01
 7.500000000000000000000000000000000000000000000000000000000000000000000000000086e-01
 8.013173181896281945892940774934203108094822908742071451762467471481287631862079e-01
 9.045084971874737120511467085914095294300772949514407155338621556763151157047318e-01
 9.251085678648070760670719614746760292353301643789564236850368282088540328017113e-01
 9.890738003669028189642833739347997662298689044313385539425888318202984165600556e-01
 9.914865498419508891409744224275993580493614375328164379986902296019539262761255e-01

In [29]:
logistic_r4_period_4_points = find_roots(x -> logistic_iterated_4(4., x), -1e-10, 1+1e-10)


Out[29]:
16-element Array{Root{Float64},1}:
 Root([-7.034439914540779e-17, 7.771423880327514e-17], :unique)
 Root([0.03376388529782207, 0.03376388529782214], :unique)     
 Root([0.043227271178699504, 0.04322727117869961], :unique)    
 Root([0.13049554138967037, 0.13049554138967054], :unique)     
 Root([0.16543469682057083, 0.165434696820571], :unique)       
 Root([0.2771308221117308, 0.27713082211173096], :unique)      
 Root([0.3454915028125262, 0.3454915028125264], :unique)       
 Root([0.45386582026834876, 0.4538658202683492], :unique)      
 Root([0.5522642316338267, 0.5522642316338269], :unique)       
 Root([0.6368314950360413, 0.6368314950360415], :unique)       
 Root([0.7499999999999999, 0.7500000000000002], :unique)       
 Root([0.8013173181896281, 0.8013173181896283], :unique)       
 Root([0.9045084971874736, 0.9045084971874738], :unique)       
 Root([0.925108567864807, 0.9251085678648071], :unique)        
 Root([0.9890738003669027, 0.9890738003669028], :unique)       
 Root([0.9914865498419508, 0.991486549841951], :unique)        

These intervals do contain the exact roots:


In [31]:
all([logistic_exact_period_4_points[i] in logistic_r4_period_4_points[i].interval for i in 1:16])


Out[31]:
true

Note that this calculation does not work using standard floating point for the calculation of the exact points -- rounding errors result in floating point numbers that are outside these guaranteed bounds.